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Abstract 

The collective modes of a rotating Bose-Einstein condensate confined in an attractive 
quadratic plus quartic trap are investigated. Assuming the presence of a large number 
of vortices we apply the diffused vorticity approach to the system. We then use the sum 
rule technique for the calculation of collective frequencies, comparing the results with the 
numerical solution of the linearized hydrodynamic equations. Numerical solutions also show 
the existence of low-frequency multipole modes which are interpreted as vortex oscillations. 

1 Introduction 

Quantized vortices are one of the most striking features of superfluids, ranging from liquid 
helium to superconductors. Bose-Einstein condensates of alkali gases have proved to be one 
of the best tools to study these fascinating quantum objects, presenting several advantages 
with respect to both helium and superconductors. On the other hand, due to the limited 
resolution imposed by current experimental techniques, direct in situ imaging of vortices 
is a very difficult task. This problem has been however overcome by using time-of-flight 
absorption images, which take advantage of the expansion of the gas after release of the 
trap. 

In this context, particularly appealing structures are given by vortex arrays, where 
singly quantized vortices typically arrange in highly regular triangular lattices, similar to 
the Abrikosov lattice of superconductors. In order to realize such configurations, large 
amounts of angular momentum have to be transferred to the gas, what can be done by 
using various experimental procedures |T|. The acquired angular velocity tends then to 
enlarge the rotating cloud, the centrifugal force giving rise to bulge effects which flatten 
the density profile towards a 2D configuration. In the presence of purely harmonic con- 
finement characterized by a frequency ojx in the plane of rotation, this phenomenon fixes 
an upper limit for the angular velocity Q. of the system, namely the frequency Q. = u±_ 
at which the quadratic effective potential given by the centrifugal force exactly equals the 
harmonic trapping term. Beyond this angular velocity the centrifugal force dominates over 
the confinement and the system is no longer bounded. 

The possibility of reaching arbitrarily high angular velocities is however provided by 
stronger than quadratic traps The introduction of a quartic term in the potential 

then opens up new regimes in the study of rotating condensates, making possible the 
realization of new equilibrium configurations with different vortex arrays. 
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In the present contribution, we will first briefly summarize the stationary solutions 
for an effectively two-dimensional Bose-Einstein condensate rotating in a harmonic plus 
quartic trap at zero temperature 0|. Then, within the Thomas- Fermi approximation 0, 
we will focus on the analysis of the most important collective excitations of the system 
|7|. Numerical solution of the hydrodynamic equations will be combined with the sum rule 
method. Vortex arrays will be treated within the so-called diffused vorticity approximation 
|S], thereby neglecting the microscopic motion of single vortices in favour of the macro- 
scopic dynamics of the system. Nevertheless, signatures of vortex modes will also be found, 
although the corresponding predicted frequencies are not expected to be very accurate. In- 
deed, a consistent calculation of such effects would require a more detailed treatment, as 
the one of Ref. 13 • 



2 Stationary configurations 

As anticipated above, in the presence of a large number of vortices, as always assumed 
in the following^ , it is possible to use the so called diffused vorticity approach, consisting 
of averaging the velocity field v over regions containing many vortex lines and assuming 
that the vorticity is spread continuously in the fluid. For example, in the case of a uniform 
vortex lattice with average vortex density n„, this corresponds to assuming rigid body 
rotation t) = f2 A r with Q, = hn^jlM, where M is the atomic mass. The relation between 
the effective angular velocity f2, oriented along the vortex line direction, and the vortex 
density n„ can be derived by imposing that the circulation of the averaged velocity field 
around a single vortex cell be equal to h/M , as for the single vortex case. More generally, the 
usual irrotationality condition V A u = of superfiuid hydrodynamics is broken in favour 
of V A v{r) — hnv{r)/M , where n„(r) is the average vortex density in the proximity of 
point r. It is then clear that the validity condition of this approach is that the average 
distance l/y'n^ between vortices be much smaller than the size of the cloud^, the concept 
of diffused vorticity being adequate to describe the dynamics only at distances larger than 
\l yfn^- If in addition ^ <^ l/y/nZ, where ^ = h/^/2gn is the healing length defined in terms 
of the interaction coupling constant g and bulk density n, one can safely use the Thomas- 
Fermi approximation. Indeed, if the size of vortex cores fixed by ^ is much smaller than 
the inter- vortex distance, one can assume a slowly varying density profile between vortices, 
consequently neglecting density gradients associated with quantum pressure effects. 

Within the presently discussed approximations, the system is then described by the 
rotational hydrodynamic equations, which in the rotating frame read 

^ + V.(™') = 0, (1) 



where n(r, t) is the spatial density, v'{r, t) — v{r ,t) — A r is the velocity in the rotating 
frame, and Vcxt is the external potential. It is trivial to check that vo = ft A r, gno = 
IJ- — Vc^t + M\fl A /2 is a stationary solution for the system, where /i, fixed by the 
normalization condition f ndr = N, is the chemical potential in the rotating frame. 

Although equilibrium configurations in Thomas-Fermi approximation can be obtained 
also for the 3D case 0, for simplicity we will consider only 2D configurations, which 
makes the analysis of collective excitations considerably easier. In fact, due to the repulsive 

^For very large rotation rates the external potential considered here is predicted to give rise to giant vortex 
states 121, where all the vorticity is confined in a single hole. We do not discuss such configurations and always 
deal with the case where singly quantized vortices are present. 

■^For an annular structure one has to compare the average vortex distance with both the annulus radius and 
the annulus width. 
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centrifugal term — |n Arp/2 in Eq. ||5J, which tends to flatten the equilibrium density no, 
this is a natural approximation for fast rotating condensates'^. Instead of using the 3D 
coupling constant gsn ~ 4Tvh?a/M, where a is the usual 3D s-wave scattering length, we 
then introduce an effective 2D coupling constant §20 ~ gao /Z, where Z is & proper length 
taking into account the extension of the real system along the rotation axis*. 
The trapping potential is given by 

huj \ / , r* 



where uj± is the harmonic oscillator frequency, dx = v h/Mui±_ is the characteristic har- 
monic oscillator length with the atomic mass M, r — y/x^ + y^ is the two-dimensional 
radial coordinate, and A is the dimensionless parameter characterizing the strength of the 
quartic term. In the following, we will use dimensionless harmonic oscillator units, where 

and d± are the units of frequency and length respectively. 

The equation for no then becomes 



1 

no = - 



" — 1 2 4 

-r — —r 



= Yg^R2-r')(r'-Rl), (4) 



2 2 

where g is the dimensionless coupling constant and 



2 - 1 , fn'^-lV 2/i 



The density is assumed to be zero where the right hand side of Eq. (|1J is negative. For 
fj. > the value of Ri becomes purely imaginary and the density vanishes at the radius 
R = R2, while for < two different radii R1.2 are present. This reflects the transition 
occurring at ^ = 0, where a hole forms in the centre of the condensate and the density 
profile assumes an annular shape. From the normalization condition one can calculate the 
transition angular velocity Qh obtaining = 1 -I- (12A^5rA''/7r)'^'^^ 0. 

For the case Q < ilh, when the hole is absent, expressing _Ri in terms oi R2 = R as 
Jii = {Q^ — 1) / \ — R? , the normalization condition gives the following third degree equation 
for R^ 

r\A\R' -in' + , (6) 

TT 

while the chemical potential becomes ^ — R?{\R' — + l)/2. In the following we will 
also need the expectation values (r^) = / nor'' dr/N = nR^{?,\R' ~ 2^' + 2)/24giV and 
(r") = TvR^iSXR' - + 5)/120gN. 

For the case > Qh, defining Rj- = R'^Ri the normalization condition simply gives 

XRt = 1^ , (7) 

TT 

and hence Jl+ — {Q,' — 1)/A and i??. — — 1)/A. The chemical potential is now fi = 
— A(_R* — R't)/8 and the previously defined expectation values become {r') = Rl/2 and 
(r*) = {5RX +R'L) /20. 



^For small angular velocities one needs the additional assumption that a strong confinement in the axial 
direction is present. 

"^For a system uniform along the 2-direction Z corresponds to the vertical size, while in the case of strong 
axial harmonic confinement one has Z = \f2iia^ , where is the oscillator length in the same direction. 
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3 Collective modes 



The collective oscillations of the system in the Thomas- Fermi approximation can be found 
by linearizing the hydrodynamic equations Q and which become 

^ <5n + V • (710 Sv) = , (8) 

^5v + g'VSn + 2nASv = . (9) 
ot 

These equations can be solved by expressing the radial and azimuthal components 5vr 
and Sv^ of the velocity field in terms of 5n and looking for solutions of the form 5n = 
5n{r)e''"^^e~"^^ , where m is the azimuthal quantum number, (j) is the azimuthal angle and 
u) is the excitation frequency in the rotating frame. This gives^ 

(a;' - 'in')SVr = i (^-Lodr + j qSu , (10) 

(lj^ - 4n^)5v^ = (^-2ndr + ^) gSn , (11) 



2 , 



2 ^^2 m~gno 2mO, d{gno) lu d ( dSn 



5n ^JL^5n+-—(rgno—]=0. (12) 

r or r or \ or ' 



At first sight, one could expect that the last equation depends both on gN and A. Actually, 
only a given combination of these parameters really matters: in particular, for given Q and 
m, the solution is uniquely fixed by the value of Qh, i-e. by the product X^gN. 

Eq. p2^ can be significantly simplified in the large Q limit 0, where useful analytical 
results can be obtained. In general, however, this equation has to be solved numerically, 
what can be achieved by direct integration with the natural initial condition® 



dSn i^{u;''-4n')+2mn\{Rl-Rl) 

-d^^^'^ - c.Ai?2(i?i-i?r ^""^^'^ ^^^^ 



and by varying u in order to obtain a well behaved solution. This procedure can be easily 
automatized by checking the validity of a condition similar to Eq. IjlH^ for the final inte- 
gration point r = for i} < Qh and r — Ri for Q > flh and essentially corresponds to 
the so called shooting method described in Ref. The code has also been checked lll| 
against the relaxation method |1()) . 

To obtain analytical results also below the large Q limit, we will rely on the sum rule 
method. To this purpose we introduce the p-energy weighted moments 

n 

relative to a generic excitation operator F = X^feLi fi^k), where Eno is the energy difference 
between the excited state \n) and the ground state |0), and a-„{F) — \{n\F\0)\^ is the 
associated strength. We also define 

m±(F) = mp{F) ± mp{F'') . (15) 

Notice that for hermitian operators F = F^ one simply has mp{F) — 2mp{F) and 
m-(F) = 0. 



^It is worth noticing that not all the solutions of Eq. 1121 correspond to physical density variations. Indeed, 
in general one has to check that the resulting eigenfunctions preserve the density normalization J &n dr = and 
that the corresponding velocity variations are finite. The latter condition, for example, leads to the exclusion of 
the solutions with lo = 2f2. 

®In fact, density boundaries are regular singular points of Eq. 1121 . For these points the second derivative 
term cancels. Notice also that one can arbitrarily fix the value of (5n(i?2), this choice being equivalent to imposing 
the amplitude of the oscillation. 
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The energy weighted moments can be used to derive rigorous upper bounds for the 
excitation frequencies of the system |12| . For example, one has the following inequality for 
the lowest energy humin excited by the operator F 



m+{F) 



(16) 



where s is positive and the equality holds whenever F excites one single mode. The explicit 
calculation of the moments ?w^^]^(_F) and m^p(_F) for p > can be carried out in terms 
of commutators between the excitation operator and the total Ifamiltonian of the system, 
evaluated on the ground state. In addition, one has the important result mlj^(_F) = — X-f(O), 
which relates the useful inverse energy weighted moment to the static limit of the dynamic 
response function xf{(^)- 

We first consider the lowest axisymmetric (m = 0) mode. One expects that this breath- 
ing oscillation is mainly excited by the monopole operator Ai = X^^Li ^i, although such 
perturbation, due to the presence of the quartic term in the potential, slightly couples also 
to higher modes^. Then, as usual, we extract the frequency of the lowest excited mode 
from the ratio between the energy weighted (mi) and inverse energy weighted (m_i) mo- 
ments. Indeed, as evident from Eq. 1141 . low order moments minimize the contributions 
coming from higher eigenfrequencies* . The mi moment for a hermitian operator can be 
expressed in terms of commutators as m\{F) = {Oj[F, [H,F]\\Q) /2. Here the many-body 
Hamiltonian H in the rotating frame, with an obvious meaning of the symbols, is given by 
H — -f/kin + ^^cxt + Hint — i^Lz, whcrc the interaction term is Hint ~ g X^i<j ^(^i ~ ^j )- For 
the monopole operator one has [M, IH,M]] = 2M and hence Tni{M) = 2N{r''). On the 
other hand, since adding a static monopole perturbation to the Hamiltonian is equivalent 
to renormalizing the trapping frequency |12| . the monopole static response can be calcu- 
lated from S{r^) = xm{0)MSlij^/2N = {d{r^) /duj'^)Suj^ (in dimensional units), where 
the derivative has to be calculated at constant angular momentum. Then, recalling that 
m_i(7M) — — xa-i(0)/2, in the Thomas- Fermi approximation one finds 



where the two expressions coincide for ft = Q^. For Q > since R-\- — (^i ~t~ -^2) — 
{Q^ — 1)/A, one finds the simple result ui = \/6QF~~2. 

Sum rule and hydrodynamic results are reported in Fig. The agreement between 
the sum rule estimate and the numerical solution of Eq. (1121 is quite good, confirming 
the hypothesis that the chosen moments are essentially saturated by the lowest m = 
mode. It is worth noticing that close to the critical angular velocity Qh the Thomas-Fermi 
approximation is expected to fail, because quantum pressure effects become important in 
the inner region of vanishing density. The sharp transition shown in the sum rule result 
and the corresponding dimple in the hydrodynamic data are in fact smoothed out by the 
full Gross-Pitaevskii solution 

We now switch to excitations of the form /(r) = j-l^le""**, which carry multipolarities 
different from zero. In particular we will concentrate on the quadrupole (m = 2) operator 
Q. The situation turns out to be much more complicated than for the monopole operator 
and one has to include in the analysis a larger number of moments, treating separately the 
two regions below and above Qh- 

Concerning the case < Qh, one can proceed exactly as in the case of purely harmonic 
trapping _12^ making a simple 2-mode assumption and solving the corresponding algebraic 
system given by the m^j^, nig , mf, and m^ moments. Indeed one expects that the low 

^In the case of purely harmonic trapping the monopole operator is instead the exact one. 
*On the other side, the coupling of the monopole operator with the next m = mode can be put in evidence 
by taking the mg/mi ratio, which is significantly higher than the chosen one. 




(17) 
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Figure 1: Lowest m = mode frequency as a function of the angular velocity for A = 0.5, 
gN = 1000 (frequencies are in harmonic oscillator units). The solid line is the sum rule estimate, 
while the triangles are obtained from the numerical solution of Eq. H12|) . The dashed vertical 
line marks the critical angular velocity fi/j = 3.2935. 



order moments of the quadrupole operator are saturated by the two lowest modes, the 
coupling with higher modes being negligible. Since m^(Q) — {Oj[Q^, Q]|0) = 0, the 2-mode 
assumption implies that the strengths of the considered m = ±2 modes are equal. A simple 
calculation then shows that the resulting frequencies are 



uj{m = ±2) = - 



+ 4- 



(18) 



while for the strengths we have am^+2{Q) = crm=-2(Q^) ~ mf /[uj{m = +2)+a>(m = —2)]. 
For the explicit calculation of the and moments in the rotating frame sum rules 
give 



m+iQ) = {0\[Q\[H,Q]]\0)^8N{r^) , 

m2{Q) = {0\[[Q\H],[H,Q]]\0) = -16N{2n{r^) 



(19) 
(20) 



where £z — —id/d<j). Since in the Thomas-Fermi diffused vorticity approach the equilibrium 
velocity is vq — ^ /\r , so that (fz) = S7(r^), one simply finds m^/m^ — —20,. In the same 
approximation, from the static quadrupole response one has instead mj^j — nR^/Sg and 
hence m+ /mt-^ = 3Ai?^ - 20,'^ + 2. Finally 



uj{m = ±2) = \/3XR^ - 02 + 2 =F ^ 



(21) 



At n = the two lowest m = ±2 modes are degenerate, but, as soon as some vorticity 
enters the system, a splitting between the modes arises^. In the rotating frame, one can 
then distinguish between a low-lying and a high-lying branch, with azimuthal quantum 
number m = +2 and m — —2 respectively. 

When O > Q,h, both the low and high-lying branch acquire an additional mode, as 
shown by the hydrodynamic numerical results reported in Fig. |21 These new modes, which 
have opposite azimuthal quantum numbers with respect to the old ones, do not have 



^Notice that at low angular velocities, when only a small number of vortices is present, the validity conditions 
of the diffused vorticity approach are not satisfied, so that Eq. 1211 is only a rough approximation. At C = 0, 
however, it gives the correct Thomas-Fermi result. 
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Figure 2: Excitation frequency as a function of the angular velocity for the main m — +2 
(upward triangles) and m = —2 (downward triangles) modes for the same parameters as in 
Fig-Gl The dashed vertical hne corresponds to = fi/i, while the solid lines correspond to the 
sum rule predictions discussed in the text. 



any counterpart for Q < Qh and arise due to the annular structure of the condensate. 
Their frequencies turn out to be very close to those of their previously discussed partners, 
becoming exactly degenerate in the large il limit [7|. Hence, since a full treatment of the 4- 
mode system would be very complicated, we base our sum rule analysis on the assumption 
that two doubly degenerate energy levels are present. In order to calculate the resulting six 
unknown quantities ojh.Li ''"h,l(Q), and an.hiQ^), one needs two more moments, namely 

{Om\H],[H,[H, Q]]]|0) = 

167V[(6n2 + l){r^) + (/) + 3A{/) - 6^(4)] , (22) 
{Q\m\H],H],[H, [H,Q]]]\0) = 
= -64iV[2t7(2f7^ + l){r^) + 2fi(/) - (60^ + + 

+3A(2n(/) - (r^l,))] , (23) 



mt(Q) 



where p 



2 



The solution of the corresponding algebraic system gives 

(24) 















fan" 


)^4< 










' rnZi 







<^h,l(Q^) 



= ± 



1 {m^ — 



f)ijJH,L + 



2 



(25) 
(26) 



where, by using the Thomas-Fermi results {l^) — Q,{r^), {p^) = Q^{r^) and mj 



Rt)/3g, one has m+/ml;^ = 4(f2^ 
1 + (3/5)A^7li/(!^2 - 1) and /r 

one can then rewrite the frequencies as 



^ en^ -2 + {6/5)x^Rt/{n^ 



n{Rl 



mJ^ /ml = 5f2 — 
- 1). Recalling that 



2 



= 30 - l-f 



3 (St - 1)' 

5 - 1 



30^^ 
5 



4 5(!^2 _ 1)2+3(^2 _ 



5 3(^2 



1)2 + {nl 



1)2 



(27) 
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Figure 3: Behaviour of the m = +2 vortex mode with the lowest number of radial nodes for the 
same parameters as in Fig.^ (a) excitation frequency w as a function of the angular velocity Q 
(the dashed vertical line corresponds to fl — ilh); (b) radial dependence of the density variation 
at n = 0.9 rih {oj = 0.396), the angular dependence being simply given by e*^"^. In the inset, the 
azimuthal component Sv^ of the velocity variation is shown for the angles (j) — 0, 7r/2, tt, Stt/2. 



The results predicted by Eqs. I|21|l and l)27|l are reported in Fig. |21 An analysis of 
the strengths given by Eqs. 12511 and 12611 shows that at SI = flh only the old modes are 
significantly excited. However, with increasing Q, while the m = +2 high-lying mode still 
has vanishing strength, the m — —2 low-lying mode becomes more and more important, 
eventually even overcoming the contribution given by the high-lying m — —2 mode |7|. It 
is also worth noticing that in the large Q limit the high-lying frequency is essentially given 

by -^J mj /to J and one has = 6S1^ — 2 as for the monopole mode. In the same limit, 

the low-lying frequency is given by o^l = (\/2/3)(r2j; — 1)/^ [7|. If one had used the same 
2- mode assumption discussed for < fih also for Q > Q^, the resulting values would have 
largely underestimated the correct frequencies. Finally, we notice that the same procedure 
could be used to extract the dipole frequencies. 

In the last part of this section we are going to discuss another class of modes found from 
the numerical solution of the hydrodynamic equations. Indeed, numerical calculations show 
the presence of very low frequency multipole modes which cannot be interpreted in terms 
of sound propagation. One can still label these excitations with the azimuthal quantum 
number to and the number of radial nodes. The frequencies of the to — +2 modes with the 
lowest number of radial nodes are plotted in Fig. OJa) as a function of the angular velocity. 
It emerges that these modes are present on both sides of the critical angular velocity 
^Ih and that their frequency is zero for f2 = 0. The rotational origin of these excitations 
together with their very low frequencies induces to identify these oscillations as multipole 
vortex modes, belonging to the family of Tkachenko modes already studied in harmonically 
trapped condensates [1411151 IT(11ll7) . In order to further investigate this hypothesis one can 
calculate the corresponding velocity variation, which indeed, at least for Q < Qh where 
the central hole is absent, resembles the typical lattice distortions of Tkachenko modes. 
This is shown in the inset of Fig.l^b), where the azimuthal component 5v^ of the velocity 
variation is plotted, the radial component Svr being practically negligible. 

Concerning the multipolarity of the modes, a simple remark is in order here. As noted 
in Ref. |18| . rotational hydrodynamic equations admit a class of zero energy solutions which 
may be interpreted as Tkachenko modes. More in detail, any zero frequency solution of 
Eqs. JHJ and Q obeys the relation 

oo 

^''^Jp^ '^sSn , (28) 
which implies ^ ■ Sv — 0. Hence, substituting into Eq. ||HJ and using dno/dcj) = 0, valid for 



8 



an axisymmetric equilibrium density profile, one finds 

|5n = 0, (29) 

so that these modes must correspond to m = 0. It follows that all the m 7^ modes, 
included Tkachenko ones, cannot have zero frequency in the diffused vorticity approach. 
However, due to the crude approximation used to treat vorticity within this method, one 
does not expect the predicted frequencies to be very accurate. 

Actually, it turns out that the numerical frequencies decrease by increasing the number 
of radial nodes. This is probably related to the fact that the nodes accumulate in the 
proximity of the cloud boundary, where the Thomas-Fermi approximation is expected to 
fail. This situation is similar to the case of the m — breathing mode at the critical angular 
velocity Qh, where the frequencies plotted in Fig.Qshow an unphysical dimple not present 
in the full Gross-Pitaevskii simulations. It is also worth noticing that the absence of such 
dimple in the positive m modes of Figs. |5] and |3^a) is indeed due to the fact that these 
oscillations are concentrated on the external boundary (see Fig.j^Jb)), so that the quality 
of the approximation used to treat the central density at f2 = f2h does not affect their 
frequency. 

According to the proposed picture, the same class of solutions of the rotational hydro- 
dynamic equations must be present also for a purely harmonic potential. This is indeed 
the case. For the m = +2 mode with the lowest number of radial nodes, where the pre- 
dicted frequency is expected to be more reliable, in the frame rotating a.t Q — 0.7 one 
finds u! = 0.204 (in units of the harmonic trapping frequency). Note that for 2D harmonic 
trapping the usual multipole modes in the rotating frame can be found analytically ac- 
cording to the formula a;(±|m|) = ^2|m| — (|m| — =p il, identical to the 3D result 
found in Ref. |18|. For the harmonic case, in addition, once the radial distance is expressed 
in units of the Thomas-Fermi radius^", Eq. H12|l depends only on m and SI, so that all the 
eigenfrequencies are independent of the interaction. The diffused vorticity estimate of the 
multipole Tkachenko frequencies, consequently, cannot properly include the compressibil- 
ity effects which have already proven to be important for the m = case Nevertheless, 
the comparison with the calculations available in the literature |19| shows that the order 
of magnitude of the predicted frequencies is correct. 



4 Conclusions 

In this paper we have presented the numerical solution of the two-dimensional linearized 
rotational hydrodynamic equations and an accurate sum rule analysis for the monopole 
and quadrupole modes of a rotating condensate in a harmonic plus quartic trap, offering 
a discussion complementary to the work contained in Ref. [J). The numerical results for 
the monopole and high-lying quadrupole modes have revealed shortcomings in using the 
Thomas-Fermi approximation rather than the full Gross-Pitaevskii solution |7| near to 
the critical angular frequency flh for hole formation. The frequency Qh has also been 
identified as the threshold angular frequency where additional modes appear, due to the 
new geometry of the system. On the other hand, sum rules have provided reliable analytical 
estimates for the frequencies and for the excitation strengths of the considered modes. 
Finally, it has been shown that the multipole vortex oscillations have a non-zero energy 
counterpart in the diffused vorticity approach, in contrast to the m = Tkachenko modes. 
The corresponding frequency estimate is however expected to be scarcely precise, due to 
the same effects which lower to zero the energy in the m = case. 

I warmly thank B. Jackson, A.L. Fetter, and S. Stringari for their precious suggestions. 

^''in physical units, the Thomas-Fermi equilibrium density for the 2D harmonic oscillator is griQ = {M/2){uj^ — 
f72)R2(l-r2/-R2). 
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